A method for atomistic spin dynamics simulations: implementation and examples 
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We present a method for performing atomistic spin dynamic simulations. A comprehensive sum- 
mary of all pertinent details for performing the simulations such as equations of motions, models 
for including temperature, methods of extracting data and numerical schemes for performing the 
simulations is given. The method can be applied in a first principles mode, where all interatomic 
exchange is calculated self-consistently, or it can be applied with frozen parameters estimated from 
experiments or calculated for a fixed spin-configuration. Areas of potential applications to different 
magnetic questions are also discussed. The method is finally applied to one situation where the 
macrospin model breaks down; magnetic switching in ultra strong fields. 



I. INTRODUCTION 



With the increasing interest in advanced magnetic ma- 
terials for data storage and processing there is an increas- 
ing need for a detailed atomistic description of magnetic 
materials. Methodological and computational schemes 
for performing atomistic magnetization dynamics have 
been presented by several groups in the pasti^^ At 
this stage however, there has not been much simulations 
on realistic systems in a materials research scope. Part 
of the reason is the computational complexity of these 
simulations. This limitation is however gradually being 
overcome by the increasing availability of computational 
power. At this stage approximate simulations of realistic 
systems are already feasible and in the future an increas- 
ing importance of atomistic modeling of magnetization 
dynamics can be expected. With recent developments in 
experimental techniques for studying magnetization dy- 
namics on short time scales and with recent findings on 
ultrafast magnetization dynamics^, there is also an in- 
creasing amount of experimental results on microscopic 
magnetization dynamics. 

The commonly used approach for studying magnetiza- 
tion dynamics, micromagnetism, provides a framework 
for understanding magnetization dynamics on length 
scales of micrometers and has with increasing compu- 
tational power become a field of large technological im- 
portance. The approach, however, suffers from a num- 
ber of limitations. It is based on the phcnomenologi- 
cal Landau-Lifshitz- Gilbert (LLG) equation where mag- 
netism is treated as a continuum vector field on a length 
scale of micrometers and where energy dissipation from 
the system is described in terms of a single ad hoc damp- 
ing parameter. This foundation limits the applicability 
and accuracy of the approach making it inadequate for 
describing various modern experiments on magnetization 
dynamics. Instead it would be desirable with an atom- 
istic approach based on the quantum description of solids, 
an approach which properly displays the connection be- 
tween the electronic structure of the material and the 
magnetization dynamics. Such an atomistic approach 
would be capable of giving a much more accurate de- 
scription of magnetization dynamics and would provide 



a framework for including a detailed description of the 
different dissipation processes involved in magnetization 
dynamics. It would provide a way of calculating magne- 
tization dynamics starting from first-principles enabling 
the study of dynamics of materials with complex chem- 
ical composition and materials with complex magnetic 
ordering such as anti-ferromagnets, spin-spirals and spin- 
glasses. 

A formal platform with which to develop an ab-initio 
spin-dynamics simulation method is naturally based on 
density functional theory, since it is known to reproduce 
both magnetic moments as well as exchange interactions 
with good accuracy. In this paper we have indeed utilized 
the efficiency of density functional theory in calculat- 
ing interatomic exchange interactions. The method pre- 
sented here is based on a Born-Oppenheimer like approx- 
imation for the spin system, where we consider the atom- 
istic spins as being slow variables, and the electronic mo- 
tion being very fast. With this adiabatic approximation 
one can separate the spin system from the electronic one, 
as shown by Antropov et al&, and hence solve the equa- 
tions of motion for the two systems separately. This is the 
approach we will adopt here but it is worth mentioning 
that alternative computational schemes for spin dynam- 
ics on an electronic level, is the time-dependent spin den- 
sity functional theory (TD-SDFT)£ or time-dependent 
current density functional theory (TD-CDFT)£. These 
approaches are promising but they are computationally 
much too time consuming for simulating larger systems. 

The scope of this article is to give a detailed presen- 
tation of a methodological and computational scheme 
for performing spin dynamic simulations on an atomistic 
scale, where most of the conceptual details were derived 
m Ref. @. The approach is hence based on an atomic 
scale description of the magnetization of a solid. Mag- 
netic properties extracted from such a description have 
long been limited to ground state properties as in density 
functional theory (DFT) or to thermal equilibrium prop- 
erties as accessed by a combination of DFT and Monte 
Carlo (MC) simulations. In Section II the adiabatic equa- 
tions of motion for the atomic spins are derived and here 
we also discuss magnetic relaxation. In Section IIIII we 
present a scheme for simulating finite temperatures in 
spin dynamics. Section llVI presents methods of extract- 
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ing and comprehending results from magnetization dy- 
namics simulations. Finally in Section IVI1 as a demon- 
stration, the method is applied to magnetic switching of 
bec Fe in ultra strong switching fields. 



II. EQUATIONS OF MOTION 
A. Slow variables 

A detailed derivation of how the dynamics of fast vari- 
ables (electrons) and slow variables (atomic spins) is sep- 
arated can, as mentioned, be found in Ref. la. Here we 
give a short description of the essential aspects of the 
dynamics of the atomic spins. 

The derivation is based on the Kohn-Sham (KS) 
Hamiltonian (Hks) of density functional theory^ which 
can be expressed in terms of a non- magnetic (Jt? nm ) and 
a magnetic part (B) as 



■J^ks = ^Sm + o ■ B. 



(1) 



The equation of motion for the slow variables, or the di- 
rections of the atomic spins, can be derived by evaluating 
the commutator between the spin operator, S, and the 
Kohn-Sham Hamiltonian 



as 

~dt 



in 



which results in 



BS 1 

— = - 7 SxB + -[S,^ nm ] 
Bt in 



(2) 



(3) 



In absence of spin-orbit coupling, S commutes with all 
terms of M^m except for the kinetic term — Yli 
(see Ref. flOl). We define the current operator as 



N 



i2m 

i 

and the spin-current operator as 
J = a ® j, 



(4) 



(5) 



where summation in Eq. ^ is performed over electrons. 
Evaluating the last term in Eq. © using the stated def- 
initions results in, 



in in \ 



^ 2r 



M]=V-J, (6) 



and the continuity equation for the spin magnetization 
within the KS framework is obtained by inserting this 
result in Eq. Q . By calculating the expectation value of 
|j| for the KS ground state we obtain, 

as, 



The second term on the left hand side is omitted for the 
applications considered in this paper. Among effects that 
arise from this term are fluctuations of the size of the 
atomic spins. For experiments when current induced ef- 
fects are important one must include this term. By using 
the atomic moment approximation (AM A), integrating 
Eq. ([7]) over atom i, we are left with a simple equation 
for the time evolution of the atomic spins, 



BS 

-^i(f) = - 7 S i (t)xB il 



(8) 



Bt 



(r, t) + V ■ J(r, t) = - 7 S(r, t) x B(r, t). (7) 



where i denotes atomic index and B^ the effective field 
which the atomic spin, Si, experiences. 



B. Parametrization 

An accurate approach for performing spin dynamics 
and for calculating effective fields acting on the atomic 
spins, is to perform a constrained DFT calculation at 
each time step using local constraining fields. This has 
been done for systems consisting of a few atoms by 
Ujfalussy et aL— where spin dynamics of a finite Co 
chain along a Pt(lll) surface step edge was simulated. 
While accurate, the approach is computationally fairly 
cumbersome and much can be gained by working with 
a parametrization of the KS Hamiltonian. Such an ap- 
proach was suggested by Fahnle et where a gradual 
trade off between accuracy and computational require- 
ments is possible. By using a spin-cluster expansion 
method, the effective field including exchange, magne- 
tocrystalline anisotropy, dipolar and external field con- 
tributions were parametrized. By increasing the number 
of parameters in the parameterization, accuracy was in- 
creased toward ab initio accuracy at the same time as 
the computational requirements increased. 

We adopt a similar approach, in the sense that the en- 
ergy of the system is parametrized and the dynamics is 
simulated for the parametrized Hamiltonian. Here, we 
present a general parametrization in terms of the atomic 
moments, mi, instead of the atomic spins, S^. For 3d- 
systems, the atomic moment is dominated by the spin 
moment contribution and the effect of spin-orbit coupling 
is small. For systems such as actinides, orbital moments 
are larger and the total atomic moment must be con- 
sidered in dynamical simulations of the magnetization. 
This would involve replacing Si for J.; in many of the 
expressions presented in this paper, but for simplicity we 
have not done so, since most of the examples presented 
here are for magnetism among 3d elements, where the or- 
bital moment can be neglected. In the KS Hamiltonian, 
we neglected dipolar interactions and spin-orbit coupling. 
Dipolar interactions are small and included separately in 
the parametrized Hamiltonian. Spin-orbit coupling gives 
rise to a magnetocrystalline anisotropy which also is in- 
cluded separately in the generalized Hamiltonian. The 
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effective field, B^, on each atom is calculated from 



B; 



dnii 



(9) 



The parametrized Hamiltonian is composed of the fol- 
lowing terms, 



(10) 



where for the first term, which represents interatomic 
exchange interactions, we use the classical Heisenberg 
Hamiltonian, 



1 \ - 



(11) 



where i and j are atomic indices, mi the classical atomic 
moment and Jij the strength of the exchange interaction. 
The second term in Eq. (JTUJ) represents the magnetocrys- 
tallinc anisotropy and can take several forms. For a uni- 
axial anisotropy we have the dominant contribution of 
the form, 



■-■"in 



K^2(m t ■ e K ) 



(12) 



where is the direction of the anisotropy axis and K 
the strength of the anisotropy field. The third term, 



(13) 



represents dipolar interactions. Here \x and v are coordi- 



nate indices and is given by, 



(14) 



where Rij is the distance between atomic moments i and 
j. Dipolar interactions are long range and important for 
the long wave length excitations. The interaction can 
be neglected in studies of short wave length excitations. 
For finite systems dipolar interactions lead to a shape 
anisotropy. For a thin film the shape anisotropy can be 
modeled by a term similar to Eq. (fT2")l . 



v/ 
Si* 



shape 



-^shapc(f*l ' ^shape) : 



(15) 



where e s h a po is the out-of-plane direction of the film, m is 
the average magnetic moment of the system and X s hapc 
is the strength of the shape anisotropy. The last term of 
Eq. dJUD, 



-B cx t • ^ 



m, 



(16) 



is the Zeeman term and describes the interaction of the 
magnetic system with an external magnetic field. 

In our approach we use the parametrized Hamiltonian, 
Eq. (dU]), combined with Eq. ([5]) which describes the 



time evolution of the magnetization for a system which 
is dominated by the spin moment. Parameters for the 
parametrized Hamiltonian are obtained by a mapping 
from a DFT ground state calculation. The most widely 
used approach is through the Liechtenstein-Katsnelson- 
Gubanov method (LKGM)^ which is based on the mag- 
netic force theorem where parameters are obtained from 
small angle perturbations from the ground state. At low 
temperatures, where the inter-atomic angles between the 
atomic spins are small, the parameters can be consid- 
ered accurate. For the paramagnetic state one may in- 
stead extract Heisenberg exchange parameters by means 
of the generalized perturbation method (GPM)^ for a 
disordered local-moment (DLM) state treated within the 
coherent potential approximation (CPA). This method 
provides a more accurate description of the high temper- 
ature region. 



C. Damping 

When the atomic spins evolve from the dynamics of 
Eq. energy and angular momentum dissipates via a 
range of mechanisms. The different mechanisms which 
lie behind this d amping have e.g. been studied in Refs. 
HIIlIlillSIllIlilSflElllallll. The effect of the different 
damping mechanisms is normally included by adding a 
phenomenological term to Eq. ([8]), which yields the LLG 
form, 

- = - 7 S i xB i + -S i x- 1 (17) 

where a is the the damping coefficient and Si is the size 
of the spin S^. For numerical reasons we use the Landau- 
Lifshitz (LL) form of damping term, and hence Eq. (I17p 
is replaced by 



dt 



-7S l x B, 



7— (Si x (Si xB,)). (18) 



III. FINITE TEMPERATURE MODELING 

Most of the systems we are interested in simulating 
with the here presented method can conceptually be un- 
derstood in terms of three thermodynamic subsystems; 
the spin system, the electronic system and the lattice 
(Fig. [T|) . The different reservoirs can be identified in mea- 
surements of specific heat. Each of these subsystems can 
be seen as reservoirs for energy and angular momentum. 

A division of the magnetic solid as such, into three 
thermodynamic reservoirs, is not free from complications, 
especially the division of the electronic system and the 
spin system which both are manifestations of the nature 
of electrons. It is important to note that the elementary 
excitations of the spin system carry an angular momen- 
tum of h. Any transfer of energy to or from the spin 
system must be accompanied by a transfer of angular 
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FIG. 1: The dynamic behavior of a magnetic solid can be 
understood in terms of three thermodynamic reservoirs and 
interactions that exchange energy between the reservoirs. In 
the figure we included approximate relaxation times within 
the reservoirs and between the reservoirs. 



momentum. The necessity of angular momentum con- 
servation is often a bottleneck of the transfer of energy 
between the subsystems. 

The processes that carry energy and angular momen- 
tum between the subsystems arc defined by the way the 
division is made. The total Hamiltonian for the mag- 
netic solid carries terms that mix the subsystems and 
these are the processes which arc responsible for the en- 
ergy and angular momentum exchange between the sub- 
systems. Relaxation rates between the reservoirs are as- 
sociated with the characteristic energies of the interac- 
tions that mediate the coupling between the reservoirs.— 
These time scales have been measured in experiments. 
The electron- lattice relaxation time, r e i, is of the order 
of picoseconds (ps). The spin-lattice relaxation time, r s i, 
is of the order of 100 ps and the spin-electron relaxation 
time, r es , has been found in recent pump-probe experi- 
ments to be of the order of 100 fs.— £ Relaxation within 
the spin system (r s ) and within the lattice (ti) are ex- 
pected to take place on times scales of the order of pi- 
coseconds whereas the electron-electron relaxation (r e ) 
takes place on a subpicosccond time scale. 

In order to describe atomistic spin dynamics at finite 
temperatures in simulations, the spin system must be 
coupled to a thermal reservoir in such a way that en- 
ergy may be transferred into and out of the system. We 
will start by showing how a single thermal reservoir can 
be coupled to the spin system and later generalize the 
discussion to several thermal reservoirs. 



A. One thermal reservoir 

For a discussion on stochastic and deterministic meth- 
ods of including temperature, see Ref. 0. One way of 
introducing a coupling to a thermal reservoir, which 
is adopted here, is through Langcvin Dynamics (LD), 
which is standard in finite temperature micromagnctic 
simulations . 25 ' 26 ! 27 ' 28 In our approach, excitations arc 
generated by performing classical rotations of single 
atomic spins in such a way that the energies of the 
atomic spins satisfy Boltzmann statistics. As a practical 
method, either Monte Carlo (MC) or LD methods may 
be used for obtaining a finite temperature equilibrium 
configuration. 

Thermal excitations are generated by adding a stochas- 
tic field, hi, to the effective field, B^, on each atom, i. 
The random field is assumed to be a Gaussian stochastic 
process with the following statistical properties, 

(&*,„(«)} = 0, (M*)Ms)) = 2DS^5 t ,S(t - a), (19) 

where (i and v are the Cartesian coordinates of the field 
and where D is the strength of the thermal fluctuations. 
The Kronecker deltas in Eq. (fT9|) state that the differ- 
ent Cartesian components of b; are unrelated and that 
the random fields acting on different magnetic moments 
i are independent. The Dirac delta states that the auto- 
correlation time of hi is much smaller than the rotational 
response of the system. 

As a technical note we mention that we have chosen to 
add the stochastic field to the effective field in both the 
precessional term and the damping term, resulting in the 
following equation: 

—i- = - 7 (S l x(B l +b l (t)))- 7 ^(S,x(S 4 x(B 4 +b i (t)))). 
at bi 

(20) 

Eq. (f20"|) is a stochastic differential equation (SDE) as 
opposed to regular ordinary differential equations (ODE) 
and require an interpretation rule^ In Appendix [S] we 
present a derivation of the amplitude of the stochastic 
field, D, required to achieve thermodynamic consistency. 

At equilibrium, MC and Spin Dynamics (SD) give 
identical results for a number of properties. MC can ac- 
tually be used as a way of benchmarking SD simulations. 
In Fig. [3 we plot the saturation magnetization for MC 
and SD for bec Fe versus temperature. Simulations are 
performed on a 20 x 20 x 20 bec system using four coor- 
dination shells in the Hcisenberg term. The Heiscnbcrg 
exchange parameters were calculated from first-principles 
theory and MC and SD arc seen to give identical results. 
In Fig. [3] we plot the energy distribution of the moments 
for MC simulations and SD simulations with two differ- 
ent damping parameters. These distributions coincide 
perfectly with the Boltzmann distribution. 

Differences between SD and MC are more subtle and 
appear first in comparisons of spin-correlation. Fig. 
illustrate the dynamic spin-correlation function S(q, w), 
which is described below in Section IIV C[ calculated 
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FIG. 2: Comparison of equilibrium magnetization versus tem- 
perature for a periodic 20 x 20 x 20 bcc Fe system for SD and 
MC. 
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FIG. 4: (Color online) The first three excitation peaks in 
S(q,u>) for a periodic, 20 x 20 x 20 bcc, Fe system at 100K. 
The size of the peaks in the SD simulation vary depending on 
the damping parameter. 
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FIG. 3: (Color online) Histogram of the energy distribution 
of the atomic spins of a 20 x 20 x 20 bcc Fe system within MC 
and SD simulations at 100K. The SD simulation is is done 
for two different damping parameters and data is obtained at 
equilibrium. 



for equilibrium states generated by MC and SD simu- 
lations with different damping parameters. For Fe realis- 
tic damping parameters are of the order 0.005-0.1^ The 
comparisons shows that the MC generated equilibrium 
has a larger amount of high energy /large momentum ex- 
citations and a lower amount of low energy excitations 
than the SD system with damping a — 0.01. By increas- 
ing the damping parameter in SD the excitation content 
is modified and for large enough damping parameter the 
number of high energy excitations exceeds that found for 
the MC equilibrium. 
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FIG. 5: (Color online) Same as in Fig.|4]but in the high energy 
region. Comparison of MC and SD with a = 0.1. 
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FIG. 6: (Color online) Same as in Fig.|4]but in the high energy 
region. Comparison of MC and SD with a = 1.0. 
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B. Several thermal reservoirs 

We have now described how the magnetic system can 
be connected to one thermal reservoir. In order to prop- 
erly represent the spin dynamics of a system, where the 
system can be decoupled into three thermal reservoirs as 
described in Section Hill we present a method of connect- 
ing the spin system to several thermodynamic reservoirs. 
The relaxation time between the electronic system and 
the lattice is of the order of picoseconds. Hence, a dis- 
tinction between the electronic reservoir and the lattice 
is only necessary when studying dynamics with resolu- 
tion higher than picoseconds. To describe the interaction 
between the electronic system, the lattice and the spin 
system, it is natural to propose a two-damping model. 
The intent is to capture the interaction between the spin 
system and the lattice with one damping parameter and 
to capture the interaction between the spin system and 
the electrons with a second damping parameter. A third 
parameter is also needed and describes the transfer of 
energy between the electrons and the lattice. Until there 
is more knowledge on how these parameters can be cal- 
culated, we use parameters obtained by fitting to pump 
probe experiments which display all these processes. 

We thus proceed by introducing two Gilbert damping 
terms which are added to the equation of motion for the 
atomic moments, 



<9S 

= - 7 S< x (B i+ bi(t))- 



-j^Si x (Si x (B< + b,-(t))) - 
- 7 ^S, x (Si x (Bi + bi(t))) 



(21) 



where a c and ot\ are the damping parameters which corre- 
spond to an energy transfer from the spin system to the 
electrons and to the lattice, respectively. Eq. ([2"T]) de- 
scribes how energy dissipates from the system through 
two channels. The temperature of the reservoirs are 
given by T e and T\. In equilibrium the temperature of 
all three thermodynamic reservoirs are the same, i.e. 
T c = T\ = T s . In our treatment the amplitude of the 
thermal fluctuations is given by, 



D = D C + D h 



where 



D , 



1 + (a e + a\) 2 7m s 



(22) 



(23) 



where x = e,l corresponds to electron or lattice effects. 
These amplitudes correspond to the equilibrium thermal 
fluctuations. Assuming a constant flux of energy from the 
reservoirs to the spin system, we use these amplitudes in 
our dynamic simulations. For practical simulations as- 
sumptions need to be made on the initial temperatures, 
T e and Tj, and the relaxation between the electrons and 
the lattice. Typical pump-probe experiments as reported 
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FIG. 7: Simulation of a pump probe experiment. Simulations 
were done on a 10 x 10 x 10 bec Fe system. The lower panel 
shows the assumed electron and lattice temperatures. The 
lattice temperature is constant at 300 K. The electronic tem- 
perature decays from an initial 500 K to 300 K with a decay 
time of 150 fs. The upper panel shows the normalized spin 
moment of the system. 



in Ref. |3l| can be simulated by assuming that the lattice 
is an infinitely large thermal reservoir with constant tem- 
perature T\. Further, we assume that the electron reser- 
voir is a thermal reservoir much smaller than the lattice, 
but much larger than the spin system with a temperature 
that evolves with time as T c (t) =T\ + T^init ■ exp(— i/ r ei) 
and where T c .init is the initial temperature. As an ap- 
plication of our two-damping model we address recent 
pump probe experiments^!, where the magnetization dy- 
namics following optical excitation of a Ni film has been 
interpreted in terms of the three thermal bath model. We 
consider a test system of bec Fe with four coordination 
shells, as described above, and are able to reproduce the 
trends found in experiments. Our simulated magnetiza- 
tion is shown in Fig. [7] In the top graph it is seen that 
the magnetization initially decreases and after some time 
(~ 0.5 ps) it stabilizes at a value ~ 70 % of the initial 
value. This behavior is in qualitative agreement with the 
measured data in Ref. K3l]. 



IV. EXTRACTING INFORMATION 

A challenge in practical simulations of magnetization 
dynamics is extracting, visualizing and comprehending 
results. A simulation of the time development of the 
magnetic moment, mi, of N atoms over M time steps, 
generates data on the form wP^tk) where j = x,y,z, 
i = [1, N] and t% = [1, £m]- For a typical simulation this 
amounts to an unmanageable amount of data which is 
difficult to store. In order to analyze and comprehend 
the meaning of the data it must be compressed into vari- 
ables that capture the state and evolution of the system. 
By doing this on the fly during the simulation, computa- 
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FIG. 8: (Color online) The trajectory of an atomic spin for a 
duration of 100 fs is shown. Simulations are performed on a 
10 x 10 x 10 bcc Fe system. On the left hand side the trajectory 
is shown for one atomic spin in a non-equilibrium system at 
OK where the atomic spin directions are completely randomly 
distributed. On the right hand side we show the trajectory 
of an atomic spin at a 300 K equilibrium. Simulations are 
performed for different sizes of times steps in the numerical 
method: 10~ 18 (black), 10~ 17 (blue) and 10~ 16 (red). 



tional time and storage requirements are greatly saved. 
Below we analyze in this way trajectories of the atomic 
spins, average moment, spin-correlations and the energy 
distributions in simulations of spin-dynamics. 



A. Trajectories 

In Fig. [5] we show the trajectories of individual atomic 
moments. The simulations are performed on a 10 x 10 x 10 
system of bcc Fe with periodic boundary conditions. The 
duration of all simulations are 100 fs and the damping is 
a = 0.1. On the left hand side we present a simulation 
at K where the initial spin distribution is random. Tra- 
jectories are presented for three different step sizes in the 
numerical scheme where Heun's scheme was used. With 
this scheme and for this particular simulation, step sizes 
as small as 1-10 attoscconds are required to produce ac- 
curate trajectories on a time scale of 100 fs. On the right 
hand side we present the trajectory of one atomic spin 
of a system in a 300 K equilibrium. We present simula- 
tions for three different step sizes. At finite temperatures 
individual trajectories do not carry much information be- 
cause of thermal fluctuations. By viewing snap shots or 
sequences of snap shots of the spin configuration over 
the entire system or parts of the system, valuable in- 
formation on correlations and domain formation can be 
visualized^ 



B. Average magnetic moment 

Averages are fundamental quantities of a magnetic sys- 
tem. With the data from a spin dynamic simulation 
averaging can be performed over space, time, different 
random number sequences in the Langcvin equations or 
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FIG. 9: (Color online) Relaxation in an easy-axis anisotropy. 
The upper panel shows the relaxation of a ferromagnet in an 
easy-axis uniaxial anisotropy and the second panel shows the 
relaxation of an anti-ferromagnet in an easy-axis the same 
anisotropy. For the ferromagnet we plot the evolution of the 
average magnetization and for the anti-ferromagnet we plot 
the evolution of the average magnetization of one sublattice. 



over different initial states. Thermal (ensemble) averages 
are often desired and can be calculated in different ways 
depending on the system and process. 

Space averaging over all atoms in the system gives the 
average magnetization. If all atoms are equivalent, such 
an average may be taken as a thermal average. Space av- 
eraging may also be performed separately over different 
sublatticcs or separately over sets of equivalent atoms in 
order to understand the the dynamical behavior of cer- 
tain parts of a system. This is useful when studying 
antifcrromagnets (AFM) or alloys. In Fig. [9] we show the 
relaxation of a ferromagnet and an antifcrromagnet in 
an uniaxial easy-axis anisotropy, following a sudden 45 
degree change with respect to the anisotropy axis. For 
the ferromagnet (upper panel) we show the evolution of 
the average magnetization whereas for the antiferromag- 
net (lower panel) we show the evolution of the average 
magnetization of a sublattice. In order to understand the 
switching behavior of an anti-ferromagnet the behavior of 
each sublattice and their mutual interaction plays an im- 
portant role. For future large scale simulations where the 
spatial variation of the magnetization over larger length 
scales is interesting, space averaging may be performed 
over several limited spatial regions of the system, produc- 
ing a more coarse grained picture of the spin dynamics of 
the system. Fig. [9] also shows that the switching of the 
AFM is faster than for the FM. 

Time averaging is useful for smearing out random fluc- 
tuations. In equilibrium, time averaging may be per- 
formed over unlimited time. When studying dynamic 
processes, time averaging must be performed over suf- 
ficiently short time intervals in comparison to the time 
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scale of the dynamic process. For systems such as spin 
glasses, where each atomic moment is unique, thermal 
averaging may be done by performing an averaging over 
time. In particular for spin glasses, which often are out 
of equilibrium, such a time averaging must be performed 
over sufficiently short time intervals. For systems with 
bond- or site-disorder, averaging can be done over differ- 
ent configurations of exchange parameters respectively 
magnetic atoms in the lattice. Dilute magnetic semicon- 
ductors are a manifestation of site-disordered systems. 
Among different classes of spin glasses there are systems 
possessing either bond- or site-disorder. Often a combi- 
nation of space and time averaging is useful. In space 
averaging the number of averaging terms is limited by 
the finite size of the system. Finite size effects them- 
selves have effects on the system which are interesting to 
study. For small systems space averaging may become 
insufficient and can be compensated by time averaging. 

Another type of averaging is averaging over identical 
simulations but with different random number sequences 
in the Langevin equations. For equilibrium simulations 
this is similar to time averaging. This type of averaging is 
however very time consuming since the same simulation 
must be repeated several times. The technique is best 
used in combination with space and time averaging. 

For some specific simulations one might also consider 
sampling over different initial configurations. The dif- 
ferent but equivalent initial configurations could be gen- 
erated with Monte Carlo or with Spin Dynamics using 
different initial configurations. 

C. Correlations between magnetic moments 

In addition to the trajectories or the absolute direc- 
tions of atomic moments, correlations or relative direc- 
tions between atomic moments, provide fundamental in- 
formation on the system (see Refs. I33ll34ll35l ). The corre- 
lation function can be defined as 

C k {y-v',t) = (m k (t)m k ,(0)) - (m k r (t))(m k r ,(0)), (24) 

where (...) denotes ensemble averages and may be per- 
formed according to the previous section. The first term 
on the right hand side is the overlap and contains infor- 
mation on the magnetic order of the system. 

In order to evaluate the spin wave excitation content of 
a system one may calculate <S(q, ui) by performing a space 
and time fourier transform of the spin-spin correlation, 

1 f + °° 
S k (q,u;) = ^yyq-(r-r) / e *»*c*(r - r', t)dt, 

v r.r' 

(25) 

where N is the number of terms in the summation. 
Fig. HHH show 5(q, u>) for an equilibrium system. For 
dynamical processes it is interesting to study how the 
spin wave content changes with time. Such a calculation 
was performed in Section IVT1 on bec Fe in an ultra strong 



switching field. For this process the time scale of the 
switching process was too fast to allow for an accurate 
S(q,cj) calculation at different points in time of the dy- 
namic process. The calculation was instead performed 
by taking snapshots of the configuration of the system 
at different points in time during the dynamic process. 
Each snap shot then serves as an initial state in a zero 
damping simulation where the dynamic spin-correlation 
is calculated. This procedure works as long as the system 
does not exhibit any strong spin- wave instabilities, which 
may change the spin- wave content at zero damping. 



D. Energy distributions 

In equilibrium the energy distribution of the spin- 
moments follows a Boltzmann distribution, as shown 
above. The energy of spin i is given by 

^ = -m 4 -B i + |m i ||B i |. (26) 

In this expression parallel coupling between moment and 
local effective field is set to zero energy. Fig. [3] shows 
a comparison between MC and SD for different values 
of the damping parameter. During dynamical processes 
the distribution changes. In analog with the spin wave 
content one may calculate the change in the energy dis- 
tribution at different points in times during a dynamic 
process. 

E. Direct visualization 

Perhaps the most natural illustration of a spin- 
dynamics simulation is to present a real time visualiza- 
tion of how the spins relax during the simulation. In 
Fig. Uni we present an example of this where a snap-shot 
of the spin-configuration of Mn doped GaAs, a diluted 
magnetic semiconductor, is shown. Only the magnetic 
Mn atoms are shown, the nonmagnetic Ga- and As-atoms 
are not visualized. The data is from a simulation where 
the Mn concentration was 5%. The temperature of the 
simulation is T = 100 K, which is below the ferromag- 
netic ordering temperature. The snap-shot illustrates 
the extent of correlation on short distances, whereas the 
value of the global magnetization is better obtained as a 
thermal average. 

V. APPLICATIONS 

An atomistic approach to spin dynamics is necessary 
for various classes of problems. One class of problems 
are systems at extreme conditions such as extreme exter- 
nal magnetic fields, where high energy short wavelength 
magnons are excited. This case will be treated in the 
following section. Another class of problems are sys- 
tems with complex magnetic ordering on an atomic scale 
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FIG. 10: (Color online) Snap-shot of a simulation of Mn 
doped GaAs (5% concentration). Only Mn atoms are shown, 
in blued color. The orange arrows indicate the size and direc- 
tion of the Mn atomic spins. 



which cannot be treated properly within micromagnetism 
such as dilute magnetic systems ) 36 ' 37 anti-ferromagnets 
and spin-spirals. A third class of problems concern sys- 
tems with complex chemical ordering or nano-structured 
materials ! 38 ' 39 In the following sections we present simu- 
lations of bcc Fe in large magnetic switching fields. We 
also illustrate different techniques for visualizing a dy- 
namic magnetization process. 



VI. MAGNETIC SWITCHING 

LLG theory relies on a macrospin approximation where 
the magnitude of the macrospin is assumed constant. It 
has previously been shown that for large anisotropies, 
the macrospin picture breaks down due to the appear- 
ance of spin-wave instabilities which alter the size of 
the macrospin. Cases where the macrospin approach 
breaks down may present interesting areas of application 
of atomistic spin dynamics. In this section we address 
magnetic switching in an external field. We show that 
the macrospin approximation remains valid up to very 
high switching field strengths. However, at extremely 
high switching fields, over 100 T, the approximation fi- 
nally breaks down. Here, we use this limit to illustrate 
the use of the atomistic spin dynamics method. By us- 
ing bulk Fe as a model system we describe the switching 
process from atomistic considerations. Furthermore we 
determine the size of the external field when LLG theory 
breaks down. 

Magnetic switching is the process of moving a system 
from one stable magnetic configuration to another and 
is fundamental for any system where a magnetic state is 
used for storing and retrieving information. The switch- 
ing process involves an excitation of the system followed 
by a relaxation into a new stable configuration. In this 
section we address magnetic switching induced by an ex- 
ternal field, in presence of which energy is transferred 
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FIG. 11: Excitation content of bcc Fe in a 300 K equilibrium 



into the magnetic system through the Zceman term. Af- 
ter this excitation, the system relaxes into a new stable 
configuration. 

Recent field pulse magnetic switching experiments ex- 
plore magnetic switching in field pulses of the order of 
35 Tj24 a factor of 1000 higher than field pulses in con- 
ventional switching experiments. As an external field is 
applied to the system, the spin- wave spectrum of the sys- 
tem is shifted by a Zeeman term which is either positive 
or negative depending on if the external field is applied 
parallel or antiparallel with the magnetization. Hence, 
the external field yields an excitation directly into the 
spin system. The system is then driven to equilibrium 
with the other thermal reservoirs by different damping 
processes. Both the orientation and magnitude of the 
macro-spin changes in a process which will be the focus 
of this section. For weak switching fields, the change 
in magnitude of the macrospin is negligible and the re- 
sults approach the LLG theory. For large fields (100 T) 
comparable to the exchange field (1000 T), the change 
in magnitude is significant and an atomistic approach is 
necessary for describing the process. 

At zero temperature our atomistic model coincides 
with the LLG picture since the system lacks high energy 
thermal excitations. At finite temperatures the system 
contains these high energy thermal excitations which al- 
ter details of the switching model. Let us now venture to 
an atomistic picture of the field induced switching process 
at finite temperatures. Since the main purpose is a quali- 
tative description of the switching process, the atomistic 
simulations are again performed on bcc Fe using four 
coordination shells in the Heisenbcrg Hamiltonian. Sim- 
ulations are performed on a system of 20 x 20 x 20 bcc 
cells. 

Consider the system in equilibrium at 300 K. Fig. [TT1 
shows the calculated 5*(q, lu) of the equilibrium state, 
which gives an image of the excitations in the system. 
Now consider what happens when a constant magnetic 
field of 1000 T is applied. We use an exaggerated large 
switching field in order to demonstrate the effect. As a 
first case consider the external field being applied parallel 
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FIG. 12: The position and magnitude of the q = 0.2e z (2-7r/a) 
peak as a 1000 T magnetic field is applied along the magneti- 
zation axis, i.e. z-axis, at t = 0. The figure shows the original 
peak at t < 0, the peak at t = when the external field is 
applied and the final peak (100 fs). The inset shows the final 
peak in blown up scale. 
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FIG. 13: (Color online) Magnetic switching process in a 
1000 T constant external field applied 135 degrees with re- 
spect to the magnetization at 300 K. The Cartesian compo- 
nents and the magnitude of the normalized (0 K) magnetic 
moment are shown. 



to the magnetization. The external field shifts the posi- 
tions of the excitations to higher energies. The process 
is illustrated in Fig. [12] for the q = 0.2e z (2ir/a) excita- 
tion peak. The figure presents the original equilibrium 
peak position for t < 0. At t = the 1000 T field is 
applied to the system resulting in a sudden shift of the 
peak. The external field adds to the exchange field which 
in turn increases the precessional torque on the atomic 
moments. This increases the frequency of precession and 
thereby the energy of the excitation. The system is now 
in a non-equilibrium state and relaxes with time through 
various dissipation processes to a new equilibrium. The 
damping torque brings the system to a new equilibrium 
with a larger saturation magnetization than the original 
equilibrium. One way to see this is that the external field 
increases the damping torque which in turn reduces the 
spread of the atomic moments and thereby increases the 
saturation magnetization. One can also regard the exter- 
nal field as a Zeeman contribution to the energy of the 
excitations. The damping term provides a path for en- 
ergy and angular momentum to leave the system. Since 
the energy of the excitations was increased by the Zee- 
man contribution, excitations need to be removed from 
the system to restore the 300 K equilibrium. Fig. [T^] 
shows how the excitation peak shrinks significantly with 
time, leaving a barely visible peak, also shown in the in- 
set. We see that in an atomistic picture, the magnitude 
of the macrospin changes in the applied magnetic field. 
This change is neglected in LLG theory where the system 
remains unchanged as a magnetic field is applied parallel 
to the magnetization. 

As a second case, let us now turn to a slightly more 
complex scenario - application of a 1000 T external field 
at an angle of 135 degrees with respect to the magne- 
tization. This scenario eventually leads to a switching 
of the direction of the macrospin. It is interesting how- 
ever to consider the path that the system takes to reach 
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FIG. 14: (Color online) Magnetic switching process in a 
1000 T constant external field applied 135 degrees with re- 
spect to the magnetization at K. The Cartesian components 
and the magnitude of the normalized (0 K) magnetic moment 
are shown. 



equilibrium in an atomistic picture. The z-axis is placed 
parallel with the external field and the initial magnetiza- 
tion lies along the (e y — e z ) direction. The switching 
processes for 300 K and K are illustrated in Fig. [TBI and 
Fig. [14] respectively. For the 300 K case, there is a small 
reduction of the magnitude of the average magnetization 
during the switching process. We also see an oscillation 
of the x- and y-components of the macrospin, signaling 
the precession of the magnetization. Let us now analyze 
the switching process in more detail. As the field is ap- 
plied to the system there is a splitting in the excitation 
spectrum. When the external field was applied parallel 
to the magnetization there was only a positive Zeeman 
contribution to the excitations. For an anti-parallel field 
we would find a negative Zeeman contribution whereas 
for a field applied in any other direction there is a nega- 
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tively shifted peak, a positively shifted peak and a peak 
at the original position. This splitting is illustrated in 
Fig. [15] for the q = 0.5e z (27r/a) peak. In this plot the 
original peak is at ~ 150 meV. The magnitude of the 
peaks with respect to each other depends only on the an- 
gle of the magnetization with respect to the external field. 
Since the angle is 135 degrees in Fig. [T5l it is the peak 
with the negative shift which dominates. The momen- 
tum dependence of the splitting of the excitation spec- 
trum is illustrated in Fig. [THJ Note that this splitting 
is a feature of the non-equilibrium system. The central 
spectrum corresponds to the equilibrium excitation spec- 
trum. Application of the 1000 T field has two main con- 
sequences. First, it excites the uniform motion of the 
magnetization or the average magnetization of the sys- 
tem. Secondly, within an atomistic model, it also splits 
the excitation peaks of the non-uniform magnons. Since 
the initial angle between magnetization and external field 
is 135 degrees, the lower branch dominates. This branch 
is lower in energy and at this instance energy starts being 
transferred from the thermal reservoir to the spin system. 
The process is illustrated in Fig. [T7] for the excitation 
peak q = 0.5e z (2-7r/a). The transfer of energy from the 
thermodynamic reservoir to the magnetic system leads 
to an increasing peak size. This initial energy transfer 
to the magnetic system leads to a reduction of the size 
of the average magnetization. However as the switching 
process continues the orientation of the average moment 
changes, reducing the angle between magnetization and 
the external field. This changes the relative strengths of 
the three peaks in the split. As the angle is reduced be- 
low 90 degrees the positive Zceman peak becomes largest 
implying an in average positive shift of the excitation. 
This reverses the energy transfer between the magnetic 
system and the thermal reservoir. Energy is now be- 
ing transported out of the magnetic system leading to 
reduced peak sizes and an increased saturation magne- 
tization. Within the atomistic model there are changes 
in the magnitude of the macrospin during the switching 
process. This does however not affect the precessional 
frequency of the macro moment which is only dependent 
on the size of the external field. 

If the initial angle between the external field and 
the magnetization is increased toward anti-parallel align- 
ment, the effect of shrinking of the magnetization during 
the switching process is enhanced. This is illustrated in 
Fig. HU where a 1000 T anti-parallel field is applied to 
the magnetization where the shrinking is total and ac- 
counts for the whole switching process. Note that at a 
certain point in the process the size of the macroscpin 
(sum of all atomic spins) is zero, after which it grows 
back to a saturation value again. For weaker external 
fields the effect is reduced. This is illustrated in Fig. [19] 
where a 100 T external field is applied anti-parallel to the 
magnetization and although the size of the applied field 
is much smaller, the magnitude of the magnetization is 
still heavily effected by the external field. 

To summarize this section, by simulations of mag- 
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FIG. 15: The splitting of the q = 0.5e z (27r/a) excitation peak 
as a 1000 T external field is applied at an angle of 135 degrees 
with respect to the magnetization. 
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FIG. 16: Splitting of the excitation spectrum during the re- 
laxation of the magnetic system in a 1000 T external field. 
The dispersion relationship for the three peaks in the upper 
left corner of Fig. 14 are shown, where circles represent the 
low energy peak, triangles the high energy peak and squares 
the middle peak. 



netization dynamics of atomic resolution, we have ex- 
plored magnetic switching in the limit of large external 
fields. A significant difference from LLG theory is seen 
for large fields. The results may aid understanding ultra- 
fast switching experiments with ultra large field pulses. 



VII. CONCLUSIONS 

A full account of the details of an atomic spin dynamics 
method has been given. A comprehensive description of 
all pertinent details for performing spin dynamics simula- 
tions, such as equations of motions, models for including 
temperature, damping mechanisms, methods of extract- 
ing data and numerical schemes for performing the sim- 
ulations, have been presented. Various ways to analyze 
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FIG. 17: Evolution of the q = 0.5e z (27r/a) excitation peak 
during the switching process. At t < 0, the peak is in its 
equilibrium position with its equilibrium magnitude. At t=0, 
a 1000 T field is applied at an angle of 135 degrees with respect 
to the magnetization. The peak is split in three. The negative 
shift is dominant in magnitude and the remaining two peaks 
are too small to be visible. Transfer of energy to the magnetic 
system leads to an increase in the peak size. As the switching 
proceeds the relative weight of the peaks change leading to 
a larger weight on the positively shifted peak. At this point, 
the energy transfer is reversed and energy is transferred out 
of the system, leading to a shrinking of all peaks. At the new 
equilibrium, the only peak remaining is the positively shifted 
peak. The magnitude of the peak is reduced compared to the 
original equilibrium peak. 
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FIG. 18: (Color online) Magnetic switching process in a 
1000 T constant external field applied anti-parallel with re- 
spect to the magnetization at 300K. 
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configuration. Various switching phenomena have been 
studied, such as the dynamics of a spin-system when 
the easy axis rapidly changes direction. Furthermore, 
we have simulated the spin dynamics of a system with 
a large applied field parallel, antiparallel and at an an- 
gle to the macro spin (the sum of all atomic spins). In 
this particular system we show that the macrospin model 
breaks down. This happens when the applied field is of 
similar size as the interatomic exchange. An experimen- 
tal realization of this is possibly best obtained for nano- 
structured magnetic multilayers, since the interatomic 
exchange interaction (which in many systems influences 
the magnetic properties heavily) typically is of the same 
size as magnetic field available in the laboratory. 
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spin dynamics simulations have been presented including 
spin-correlations. The developed method can be applied 
in a first principles mode, where all interatomic exchange 
is calculated self-consistently, or it can be applied with 
frozen parameters estimated from experiments or calcu- 
lated for a fixed spin-configuration. 

Applications of the method have been made to sev- 
eral systems, with interatomic exchange calculated from 
first principles, primarily for bec Fc from a frozen spin- 



APPENDIX A: LANGEVIN SPIN DYNAMICS 

The Fokker-Planck equation describes the time evolu- 
tion of a non-equilibrium probability distribution. The 
Fokker-Planck equation corresponding to the stochas- 
tic Landau Lifshitz Gilbert (LLG) equation has been 
derived 4^ 

The general form of the Langevin equations can be 
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written as 



dy 
dt 



Ai{y,t) + Y,B ik (y,t)L k (t), 



(Al) 



This condition on the Fokkcr-Planck equation is not con- 
sistent with Ito calculus. For Stratanovich calculus, we 
can find a condition on r that makes the equations fulfill 
the equilibrium requirement. First note that 



where 

< L k (t) >= 0, < L k {t)Li{s) >= 2D5 u S{t - s). (A2) 

The stochastic LLG equation can be written in the gen- 
eral form of a Langevin-equation by identifying the coef- 
ficients 



A: = 7l m x B cff m x (m x B off )] , 

m i 



(A3) 



B ik = 7^ tijkrrij + ^-{m 2 8 lk ~ mim k )}. (A4) 



The time evolution of the general form, using 
Stratanovich calculus is given by 



d 



E^^E^)P] (A5) 



B 



off 



dJf? 
dm 



Using Eq. (|A8|) we can write 



dPo -m P 
am 

Hence, the first term of Eq. (|A6|) . 



dm 



(7m x B cff )P 



vanishes and the remainder can be written as 



dP d , . a , . 

-57 = ~ q— I -7— m x (m x B off 



2r w 



(A9) 



(A10) 



(All) 



m x (m x B eff )]P}, (A12) 



We then arrive at the following Fokker-Plank equation for 
the time evolution of the probability distribution, P(m), 
of the atomic moments, 

<9P 9 .. a . _ . 

— = - — { 7m x B cf f- 7— m x (m x B off ) + 
at am m 

1 d 
+— mx (mx — )]P},(A6) 
It am 



where we have defined t as, 
1 



2£>7 ii (l + a 2 ) 



(A7) 



The Fokker-Planck equation associated with the stochas- 
tic LLG equation must satisfy the correct thermal- 
equilibrium properties. In thermal equilibrium P(m) 
must have the form of the Boltzmann distribution, 



From this we see that a requirement for a stationary so- 
lution, or dP/dt = 0, is that ja/m = (3/2tm- Hence, 
resulting in 



TN 



1 



(A13) 



In Eq. (|A7|) we sec that the temperature determines the 
amplitude of the random field, b^, and this is how tem- 
perature enters our simulations. The amplitudes are fi- 
nally given by, 



D 



a ksT 
1 + a 2 7m ' 



(A14) 



P (m) oc cxp[-/3jr(m)]. 



(A8) 
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